#Load necessary packages and data
packages <- c("ggplot2", "dplyr", "tidyr", "Hmisc", "skimr", "car", "nortest", "lmtest", "texreg", "writexl")
purrr::walk(packages, library, character.only = TRUE)

getwd()
data <- read.csv("Supplementary Data 1 - database.csv")

str(data)

#Exclude data with missing Severity
data <- data |>
  filter(!is.na(Severity))

#Descriptive analysis
skimr::skim(data)

skim_geral <- skimr::skim(data)
skim_geral <- skim_geral |>
  mutate(across(where(is.numeric), ~ format(., scientific = FALSE)))
skim_geral_df <- as.data.frame(skim_geral)
write_xlsx(skim_geral_df, "descriptive.xlsx")


dummies <- c("PhD", "Bachelor", "BachArea", "Professor", "Papers")
contagem_valores <- sapply(data[dummies], function(x) table(factor(x, levels = c(0, 1))))
print(contagem_valores)


#Sector Ajustments
data <- data |>
  mutate(
    Agricultura = `SetorAgricultura..Açúcar..Álcool.e.Cana.`,
    Alimentos = SetorAlimentos + `SetorEmp..Adm..Part....Alimentos`,
    Bebidas_Fumo = `SetorBebidas.e.Fumo`,
    Comércio = `SetorComércio..Atacado.e.Varejo.`,
    Comunicação_Informática = `SetorComunicação.e.Informática` + `SetorEmp..Adm..Part....Comunicação.e.Informática`,
    Const_Civil = `SetorConstrução.Civil..Mat..Constr..e.Decoração` + `SetorEmp..Adm..Part....Const..Civil..Mat..Const..e.Decoração`,
    Educação = SetorEducação + `SetorEmp..Adm..Part....Educação`,
    Energia_Elétrica = `SetorEmp..Adm..Part....Energia.Elétrica` + `SetorEnergia.Elétrica`,
    Total_Farmacêutico_Higiene = `SetorFarmacêutico.e.Higiene`,
    Máqs_Equip_Veíc_Peças = `SetorMáquinas..Equipamentos..Veículos.e.Peças`,
    Metalurgia_Siderurgia = `SetorMetalurgia.e.Siderurgia` + `SetorEmp..Adm..Part....Metalurgia.e.Siderurgia`,
    Outros = SetorOutros,
    Petróleo_Gás = `SetorEmp..Adm..Part....Petróleo.e.Gás` + `SetorPetróleo.e.Gás`,
    Saneamento_Serv_Água_Gás = `SetorSaneamento..Serv..Água.e.Gás`,
    Serviços_Médicos = `SetorServiços.Médicos`,
    Serviços_Transporte_Logística = `SetorEmp..Adm..Part....Serviços.Transporte.e.Logística` + `SetorServiços.Transporte.e.Logística`,
    Telecomunicações = SetorTelecomunicações,
    Têxtil_Vestuário = `SetorTêxtil.e.Vestuário`
  ) |>
  select(
    -`SetorAgricultura..Açúcar..Álcool.e.Cana.`,
    -SetorAlimentos,
    -`SetorEmp..Adm..Part....Alimentos`,
    -`SetorBebidas.e.Fumo`,
    -`SetorComércio..Atacado.e.Varejo.`,
    -`SetorComunicação.e.Informática`,
    -`SetorEmp..Adm..Part....Comunicação.e.Informática`,
    -`SetorConstrução.Civil..Mat..Constr..e.Decoração`,
    -`SetorEmp..Adm..Part....Const..Civil..Mat..Const..e.Decoração`,
    -SetorEducação,
    -`SetorEmp..Adm..Part....Educação`,
    -`SetorEmp..Adm..Part....Energia.Elétrica`,
    -`SetorEnergia.Elétrica`,
    -`SetorFarmacêutico.e.Higiene`,
    -`SetorMáquinas..Equipamentos..Veículos.e.Peças`,
    -`SetorMetalurgia.e.Siderurgia`,
    -`SetorEmp..Adm..Part....Metalurgia.e.Siderurgia`,
    -SetorOutros,
    -`SetorPetróleo.e.Gás`,
    -`SetorEmp..Adm..Part....Petróleo.e.Gás`,
    -`SetorSaneamento..Serv..Água.e.Gás`,
    -`SetorServiços.Médicos`,
    -`SetorEmp..Adm..Part....Serviços.Transporte.e.Logística`,
    -`SetorServiços.Transporte.e.Logística`,
    -SetorTelecomunicações,
    -`SetorTêxtil.e.Vestuário`
  )

setores <- c("Agricultura", "Alimentos", "Bebidas_Fumo", "Comércio", 
             "Comunicação_Informática", "Const_Civil", "Educação", 
             "Energia_Elétrica", "Total_Farmacêutico_Higiene", 
             "Máqs_Equip_Veíc_Peças", "Metalurgia_Siderurgia", "Outros", 
             "Petróleo_Gás", "Saneamento_Serv_Água_Gás", "Serviços_Médicos", 
             "Serviços_Transporte_Logística", "Telecomunicações", 
             "Têxtil_Vestuário")


tabela_setores <- data |>
  select(all_of(setores)) |>
  summarise(across(everything(), ~ sum(. == 1))) |>
  pivot_longer(cols = everything(), names_to = "Setor", values_to = "Frequência") |>
  arrange(desc(Frequência))

print(tabela_setores)

#Main Analysis - Severity
model_0 <- lm(Severity ~ idade + log_ativo_total + alavancagem + Segmento_GCBásico + Segmento_GCNovo.Mercado + Agricultura + Alimentos + Bebidas_Fumo + Comércio + Comunicação_Informática + Const_Civil + Educação + Energia_Elétrica + Total_Farmacêutico_Higiene + Máqs_Equip_Veíc_Peças + Metalurgia_Siderurgia + Petróleo_Gás + Saneamento_Serv_Água_Gás + Serviços_Médicos + Serviços_Transporte_Logística + Telecomunicações,data = data)
summary(model_0)

model_1 <- lm(Severity ~ PhD_holder_board + Bachelor_board + Area_board + Professor_board + Papers_published_board + Bachelor + idade + log_ativo_total + alavancagem + Segmento_GCBásico + Segmento_GCNovo.Mercado + Agricultura + Alimentos + Bebidas_Fumo + Comércio + Comunicação_Informática + Const_Civil + Educação + Energia_Elétrica + Total_Farmacêutico_Higiene + Máqs_Equip_Veíc_Peças + Metalurgia_Siderurgia + Petróleo_Gás + Saneamento_Serv_Água_Gás + Serviços_Médicos + Serviços_Transporte_Logística + Telecomunicações, data = data)
summary(model_1)

model_2 <- lm(Severity ~ Intensity * PhD_holder_board +  Intensity * Bachelor_board + Intensity * Area_board + Intensity * Professor_board + Intensity * Papers_published_board +
                Intensity * Bachelor +
                idade + log_ativo_total + alavancagem + 
                Segmento_GCBásico + Segmento_GCNovo.Mercado+ 
                Agricultura + Alimentos + Bebidas_Fumo + Comércio + Comunicação_Informática + Const_Civil + Educação + Energia_Elétrica + Total_Farmacêutico_Higiene + Máqs_Equip_Veíc_Peças + Metalurgia_Siderurgia +  Petróleo_Gás + Saneamento_Serv_Água_Gás + Serviços_Médicos + Serviços_Transporte_Logística + Telecomunicações,
              data = data)
summary(model_2)

shapiro_francia <- sf.test(residuals(model_2))
print(shapiro_francia)

bp_test <- bptest(model_2)
print(bp_test)

dw_test <- dwtest(model_2)
print(dw_test)


model_list <- list(model_0, model_1, model_2)
screenreg(model_list, 
          custom.model.names = c("Model 0", "Model 1", "Model 2"),
          stars = c(0.001, 0.01, 0.05, 0.1), 
          single.row = FALSE,
          include.rsquared = TRUE, 
          digits = 3, 
          decimal.mark = ",")  

htmlreg(model_list, 
        file = "tabela_severidade.html",  
        custom.model.names = c("Model 0", "Model 1", "Model 2"),
        stars = c(0.001, 0.01, 0.05, 0.1), 
        single.row = FALSE,
        include.rsquared = TRUE,  
        digits = 3, 
        decimal.mark = ",")

#Cox Regression - Time to Recovery
library("survival")
library("survminer")
quantile(data$Time_Recovery, probs = c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE)

data$Event_Indicator <- ifelse(is.na(data$Time_Recovery) | data$Time_Recovery > 55, 0, 1)

res.cox_0 <- coxph(Surv(Time_Recovery, Event_Indicator) ~ idade + log_ativo_total + alavancagem + Segmento_GCNovo.Mercado,data = data)
summary(res.cox_0)

res.cox_1 <- coxph(Surv(Time_Recovery, Event_Indicator) ~ PhD_holder_board + Bachelor_board + Professor_board + Bachelor + idade + log_ativo_total + alavancagem + Segmento_GCNovo.Mercado, data = data)
summary(res.cox_1)

res.cox_2 <- coxph(Surv(Time_Recovery, Event_Indicator) ~ Intensity * PhD_holder_board +  Intensity * Bachelor_board + Intensity * Professor_board +
                     idade + log_ativo_total + alavancagem +
                     Segmento_GCNovo.Mercado,
                   data = data)
summary(res.cox_2)

model_list <- list(res.cox_0, res.cox_1, res.cox_2)
screenreg(model_list, 
          custom.model.names = c("Model 0", "Model 1", "Model 2"),
          stars = c(0.001, 0.01, 0.05, 0.1), 
          single.row = FALSE,
          include.rsquared = TRUE, 
          digits = 3, 
          decimal.mark = ",")  
htmlreg(model_list, 
        file = "tabela_temporecuperação.html",  
        custom.model.names = c("Model 0", "Model 1", "Model 2"),
        stars = c(0.001, 0.01, 0.05, 0.1), 
        single.row = FALSE,
        include.rsquared = TRUE,  
        digits = 3, 
        decimal.mark = ",")

#Robustnes Check
library(quantreg)

model_q10 <- rq(Severity ~ Intensity * PhD_holder_board +  Intensity * Bachelor_board + Intensity * Area_board + 
                  Intensity * Professor_board + Intensity * Papers_published_board +
                  Intensity * Bachelor +
                  idade + log_ativo_total + alavancagem + 
                  Segmento_GCBásico + Segmento_GCNovo.Mercado + 
                  Agricultura + Alimentos + Bebidas_Fumo + Comércio + Comunicação_Informática + 
                  Const_Civil + Educação + Energia_Elétrica + Total_Farmacêutico_Higiene + 
                  Máqs_Equip_Veíc_Peças + Metalurgia_Siderurgia + Petróleo_Gás + 
                  Saneamento_Serv_Água_Gás + Serviços_Médicos + Serviços_Transporte_Logística + 
                  Telecomunicações,
                data = data,
                tau = 0.1) 
summary(model_q10)


model_q20 <- rq(Severity ~ Intensity * PhD_holder_board +  Intensity * Bachelor_board + Intensity * Area_board + 
                  Intensity * Professor_board + Intensity * Papers_published_board +
                  Intensity * Bachelor +
                  idade + log_ativo_total + alavancagem + 
                  Segmento_GCBásico + Segmento_GCNovo.Mercado + 
                  Agricultura + Alimentos + Bebidas_Fumo + Comércio + Comunicação_Informática + 
                  Const_Civil + Educação + Energia_Elétrica + Total_Farmacêutico_Higiene + 
                  Máqs_Equip_Veíc_Peças + Metalurgia_Siderurgia + Petróleo_Gás + 
                  Saneamento_Serv_Água_Gás + Serviços_Médicos + Serviços_Transporte_Logística + 
                  Telecomunicações,
                data = data,
                tau = 0.2) 
summary(model_q20)

model_q30 <- rq(Severity ~ Intensity * PhD_holder_board +  Intensity * Bachelor_board + Intensity * Area_board + 
                  Intensity * Professor_board + Intensity * Papers_published_board +
                  Intensity * Bachelor +
                  idade + log_ativo_total + alavancagem + 
                  Segmento_GCBásico + Segmento_GCNovo.Mercado + 
                  Agricultura + Alimentos + Bebidas_Fumo + Comércio + Comunicação_Informática + 
                  Const_Civil + Educação + Energia_Elétrica + Total_Farmacêutico_Higiene + 
                  Máqs_Equip_Veíc_Peças + Metalurgia_Siderurgia + Petróleo_Gás + 
                  Saneamento_Serv_Água_Gás + Serviços_Médicos + Serviços_Transporte_Logística + 
                  Telecomunicações,
                data = data,
                tau = 0.3) 
summary(model_q30)



model_q40 <- rq(Severity ~ Intensity * PhD_holder_board +  Intensity * Bachelor_board + Intensity * Area_board + 
                  Intensity * Professor_board + Intensity * Papers_published_board +
                  Intensity * Bachelor +
                  idade + log_ativo_total + alavancagem + 
                  Segmento_GCBásico + Segmento_GCNovo.Mercado + 
                  Agricultura + Alimentos + Bebidas_Fumo + Comércio + Comunicação_Informática + 
                  Const_Civil + Educação + Energia_Elétrica + Total_Farmacêutico_Higiene + 
                  Máqs_Equip_Veíc_Peças + Metalurgia_Siderurgia + Petróleo_Gás + 
                  Saneamento_Serv_Água_Gás + Serviços_Médicos + Serviços_Transporte_Logística + 
                  Telecomunicações,
                data = data,
                tau = 0.4) 
summary(model_q40)



model_q50 <- rq(Severity ~ Intensity * PhD_holder_board +  Intensity * Bachelor_board + Intensity * Area_board + 
                  Intensity * Professor_board + Intensity * Papers_published_board +
                  Intensity * Bachelor +
                  idade + log_ativo_total + alavancagem + 
                  Segmento_GCBásico + Segmento_GCNovo.Mercado + 
                  Agricultura + Alimentos + Bebidas_Fumo + Comércio + Comunicação_Informática + 
                  Const_Civil + Educação + Energia_Elétrica + Total_Farmacêutico_Higiene + 
                  Máqs_Equip_Veíc_Peças + Metalurgia_Siderurgia + Petróleo_Gás + 
                  Saneamento_Serv_Água_Gás + Serviços_Médicos + Serviços_Transporte_Logística + 
                  Telecomunicações,
                data = data,
                tau = 0.5) 
summary(model_q50)

model_q60 <- rq(Severity ~ Intensity * PhD_holder_board +  Intensity * Bachelor_board + Intensity * Area_board + 
                  Intensity * Professor_board + Intensity * Papers_published_board +
                  Intensity * Bachelor +
                  idade + log_ativo_total + alavancagem + 
                  Segmento_GCBásico + Segmento_GCNovo.Mercado + 
                  Agricultura + Alimentos + Bebidas_Fumo + Comércio + Comunicação_Informática + 
                  Const_Civil + Educação + Energia_Elétrica + Total_Farmacêutico_Higiene + 
                  Máqs_Equip_Veíc_Peças + Metalurgia_Siderurgia + Petróleo_Gás + 
                  Saneamento_Serv_Água_Gás + Serviços_Médicos + Serviços_Transporte_Logística + 
                  Telecomunicações,
                data = data,
                tau = 0.6) 
summary(model_q60)

model_q70 <- rq(Severity ~ Intensity * PhD_holder_board +  Intensity * Bachelor_board + Intensity * Area_board + 
                  Intensity * Professor_board + Intensity * Papers_published_board +
                  Intensity * Bachelor +
                  idade + log_ativo_total + alavancagem + 
                  Segmento_GCBásico + Segmento_GCNovo.Mercado + 
                  Agricultura + Alimentos + Bebidas_Fumo + Comércio + Comunicação_Informática + 
                  Const_Civil + Educação + Energia_Elétrica + Total_Farmacêutico_Higiene + 
                  Máqs_Equip_Veíc_Peças + Metalurgia_Siderurgia + Petróleo_Gás + 
                  Saneamento_Serv_Água_Gás + Serviços_Médicos + Serviços_Transporte_Logística + 
                  Telecomunicações,
                data = data,
                tau = 0.7) 
summary(model_q70)


model_q80 <- rq(Severity ~ Intensity * PhD_holder_board +  Intensity * Bachelor_board + Intensity * Area_board + 
                  Intensity * Professor_board + Intensity * Papers_published_board +
                  Intensity * Bachelor +
                  idade + log_ativo_total + alavancagem + 
                  Segmento_GCBásico + Segmento_GCNovo.Mercado + 
                  Agricultura + Alimentos + Bebidas_Fumo + Comércio + Comunicação_Informática + 
                  Const_Civil + Educação + Energia_Elétrica + Total_Farmacêutico_Higiene + 
                  Máqs_Equip_Veíc_Peças + Metalurgia_Siderurgia + Petróleo_Gás + 
                  Saneamento_Serv_Água_Gás + Serviços_Médicos + Serviços_Transporte_Logística + 
                  Telecomunicações,
                data = data,
                tau = 0.8) 
summary(model_q80)

model_q90 <- rq(Severity ~ Intensity * PhD_holder_board +  Intensity * Bachelor_board + Intensity * Area_board + 
                  Intensity * Professor_board + Intensity * Papers_published_board +
                  Intensity * Bachelor +
                  idade + log_ativo_total + alavancagem + 
                  Segmento_GCBásico + Segmento_GCNovo.Mercado + 
                  Agricultura + Alimentos + Bebidas_Fumo + Comércio + Comunicação_Informática + 
                  Const_Civil + Educação + Energia_Elétrica + Total_Farmacêutico_Higiene + 
                  Máqs_Equip_Veíc_Peças + Metalurgia_Siderurgia + Petróleo_Gás + 
                  Saneamento_Serv_Água_Gás + Serviços_Médicos + Serviços_Transporte_Logística + 
                  Telecomunicações,
                data = data,
                tau = 0.9) 
summary(model_q90)


boot10 <- summary(model_q10, se = "boot", R = 1000) #bootstrap
boot20 <- summary(model_q20, se = "boot", R = 1000)
boot30 <- summary(model_q30, se = "boot", R = 1000)
boot40 <- summary(model_q40, se = "boot", R = 1000)
boot50 <- summary(model_q50, se = "boot", R = 1000)
boot60 <- summary(model_q60, se = "boot", R = 1000)
boot70 <- summary(model_q70, se = "boot", R = 1000)
boot80 <- summary(model_q80, se = "boot", R = 1000)
boot90 <- summary(model_q90, se = "boot", R = 1000)


extract_summary_rq_ready <- function(summary_obj, tau_value) {
  if ("coefficients" %in% names(summary_obj)) {
    coef_df <- as.data.frame(summary_obj$coefficients)
    coef_df$tau <- tau_value
    return(coef_df)
  } else {
    stop("Os coeficientes não foram encontrados no objeto summary.rq.")
  }
}

summaries <- list(boot10, boot20, boot30, boot40, boot50, 
                  boot60, boot70, boot80, boot90)

tau_values <- c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9)

for (i in 1:length(summaries)) {
  boot_df <- extract_summary_rq_ready(summaries[[i]], tau_values[i])
  write.csv(boot_df, file = paste0("resultados_bootstrap_tau_", tau_values[i], ".csv"), row.names = FALSE)
}


